///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/ham/ham.cpp
///
/// \brief   some implementation for the hamiltonian class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<formic/ham/ham.h>
#include<formic/archive.h>
#include<formic/input/base.h>
#include<formic/exception.h>
#include<formic/mpi/interface.h>
#include<formic/psi3/interface.h>
#include<formic/lapack/interface.h>
//#include<formic/molpro/interface.h>
#include<formic/timing/timing.h>
#include<formic/output/output.h>
#include<formic/command/command.h>
#include<formic/all_data.h>
#include<formic/ham/integral_transform.h>
#include<formic/numeric/matrix.h>
#include<formic/ham/er_local.h>
#include<formic/ham/lowdin_local.h>
#include<formic/ham/freeze_core.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   initializes an empty hamiltonian with the appropriate number of orbitals
///
/// \param[in]   norb       the number of basis orbitals
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> formic::Hamiltonian<S>::Hamiltonian(const int norb)
  : _norb(norb),
    _e_cor_nuc(0.0),
    _orb_coeff(norb, norb, formic::zero(S())),
    _ovi(norb*norb, formic::zero(S())),
    _oei(norb*norb, formic::zero(S())),
    _mei(norb*norb, formic::zero(S())),
    _tei(formic::tei_compact_array_length(norb), formic::zero(S())),
    _n_cholesky(0),
    _chol_thresh(0.0),
    _cholesky_vecs(),
    _cholesky_vecs_t(),
    _dense_tei()
{
}

template formic::Hamiltonian<double>::Hamiltonian(const int norb);
template formic::Hamiltonian<std::complex<double> >::Hamiltonian(const int norb);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   sets the values of the one- and two-electron integrals
///
/// \param[in]   new_oei    one electron integrals, size _norb*_norb
/// \param[in]   new_tei    two electron integrals, dense storage, size _norb*_norb*_norb*_norb
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void formic::Hamiltonian<S>::set_integrals(const S * const new_oei, const S * const new_tei) {

  // root process sets the integrals
  if ( formic::mpi::rank() == 0 ) {
    std::copy(new_oei, new_oei + _norb*_norb, _oei.begin());
    formic::tei_dense_to_compact(_norb, new_tei, &_tei.at(0));
  }

  // broadcast integrals to all processes
  formic::mpi::bcast(_oei);
  formic::mpi::bcast(_tei);

  // clear dense tei array (it will be rebuilt the next time it is needed)
  _dense_tei.clear();

  // re-initialize the modified one electron integrals
  this->initialize_mei();

  //this->compute_cholesky();

}

template void formic::Hamiltonian<double>::set_integrals(const double * const new_oei, const double * const new_tei);
template void formic::Hamiltonian<std::complex<double> >::set_integrals(const std::complex<double>  * const new_oei, const std::complex<double>  * const new_tei);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   transforms the one and two electron integrals based on the supplied orbital rotation
///
/// \param[in]   orb_rot      matrix of orbital transformation, size _norb*_norb, with the
///                           pth coefficient of the ith new orbital stored in orb_rot[i+p*_norb]
/// \param[in]   from_curr    whether to transform starting from the current or initial orbitals
/// \param[in]   prihnt_only  if true, will print what the new orbitals would be without actually
///                           making any changes to the Hamiltonian
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::Hamiltonian<S>::transform_orbitals(const S * const orb_rot, const bool from_curr, const bool print_only) {

  // start with either the current orbitals or the initial atomic orbitals
  const formic::ConstMatrix<S,int> start = ( from_curr ? _orb_coeff.clone() : formic::identity_matrix<S,int>(_norb) );

  formic::Matrix<S,int> new_coeff;
  if ( formic::mpi::rank() == 0 ) {

    // compute new orbital coefficients as the transformation of the current orbitals
    new_coeff = formic::Matrix<S,int>(_norb, _norb, formic::zero(S()));
    formic::xgemm('N', 'N', _norb, _norb, _norb, formic::unity(S()), orb_rot, _norb, &start.at(0,0), _norb, formic::zero(S()), &new_coeff.at(0,0), _norb);

    // print new orbitals
    if ( print_only )
      formic::of << boost::format("orbital basis:") << std::endl;
    else
      formic::of << boost::format("transforming Hamiltonians integrals into the following orbital basis:") << std::endl;
    for (int p = 0; p < _norb; p++) {
      for (int i = 0; i < _norb; i++)
        formic::of << boost::format(" %14.8f") % new_coeff.at(i,p);
      formic::of << std::endl;
    }
    formic::of << std::endl;

  }

  // if only printing what the new orbitals would be, stop here
  if ( print_only )
    return;

  // read initial integrals from disk
  this->from_archive("fqmc_init_hamiltonian.bin");

//  // print initial integrals
//  if ( formic::mpi::rank() == 0 ) {
//    this->print_oei(_norb, _oei, "initial one electron integrals:");
//    this->print_tei_diag(_norb, _tei, "[pp|qq] elements of initial two electron integrals:");
//  }

  // compute transformation on root process
  if ( formic::mpi::rank() == 0 ) {

    // set orbital coefficients
    _orb_coeff <<= new_coeff;

    // get space to hold the transformed integrals
    std::vector<S> new_oei(_norb*_norb, formic::zero(S()));
    std::vector<S> new_tei(_norb*_norb*_norb*_norb, formic::zero(S()));

    // perform the transformation
    formic::transform_dense_oei_tei(_norb, &_oei.at(0), this->dense_tei_ptr(), &_orb_coeff.at(0,0), &new_oei.at(0), &new_tei.at(0));

    // empty the dense two electron integrals (they will be rebuilt the next time they are needed)
    _dense_tei.clear();

    // update the one electron integrals
    _oei = new_oei;

    // update the 8-fold symmetry storage of the two electron integrals
    formic::tei_dense_to_compact(_norb, &new_tei.at(0), &_tei.at(0));
    //for (int c = 0; c < _norb; c++)
    //for (int d = 0; d <= c   ; d++)
    //for (int a = 0; a < _norb; a++)
    //for (int b = 0; b <= a   ; b++)
    //  this->get_tei(a,b,c,d) = new_tei[a+_norb*(b+_norb*(c+_norb*d))];

    // transform the overlap matrix
    std::vector<S> old_ovi(_ovi);
    formic::transform_oei(_norb, &_orb_coeff.at(0,0), &old_ovi.at(0), &_ovi.at(0));

    // re-initialize the modified one electron integrals
    this->initialize_mei();

    // recompute the cholesky decomposition
    this->compute_cholesky();

  }

  // empty the dense two electron integrals (they will be rebuilt the next time they are needed)
  _dense_tei.clear();

  // root writes the hamiltonian to disk
  if ( formic::mpi::rank() == 0 )
    this->to_archive("fqmc_trns_hamiltonian.bin");

  // read in and broadcast the transformed hamiltonian to all processes
  this->from_archive("fqmc_trns_hamiltonian.bin");

  // print transformed integrals
  if ( formic::mpi::rank() == 0 ) {
    this->print_oei(_norb, _oei, "transformed one electron integrals:");
    this->print_oei(_norb, _ovi, "transformed overlap integrals:");
    this->print_tei_diag(_norb, _tei, "[pp|qq] elements of transformed two electron integrals:");
  }

}

template void formic::Hamiltonian<double>::transform_orbitals(const double * const, const bool, const bool);
template void formic::Hamiltonian<std::complex<double> >::transform_orbitals(const std::complex<double> * const, const bool, const bool);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   reads the overlap matrix and one and two electron integrals from a text file
///
/// \param[in]   userinp  user input options
/// \param[in]   filename the file from which to read the integrals
/// \param[out]  enuc     the nuclear repulsion energy
/// \param[out]  ovi      the atomic orbital overlap integrals, size = no*no
/// \param[out]  oei      the one electron integrals, size = no*no
/// \param[out]  tei      the two electron integrals, size = no*(no+1) * ( no*(no+1) + 2 ) / 8
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void formic::get_text_ints(const formic::InputBase & userinp,
                                             const std::string & filename,
                                             double & e_core_nuc,
                                             S * const ovi,
                                             S * const oei,
                                             S * const tei) {

  // open the file
  std::ifstream ifs(filename.c_str());
  if (!ifs.is_open())
    throw formic::Exception("failed to open file \"%s\" in formic::get_text_ints") % filename;

  // read the number of orbitals
  int no = -1;
  ifs >> no;
  if ( no != userinp.get<int>("nao") )
    throw formic::Exception("wrong number of orbitals \"%i\" in formic::get_text_ints") % no;

//  // read the number of irreps
//  int nirr = -1;;
//  ifs >> nirr;
//
//  // read the orbital of each irrep
//  std::vector<int> opirr(nirr, 0);
//  std::vector<int> orb_irrep(no, 0);
//  for (int i = 0; i < no; i++) {
//    ifs >> orb_irrep.at(i);
//    opirr.at(orb_irrep.at(i)) += 1;
//  }

  // read the core and nuclear energy
  ifs >> e_core_nuc;

  // read the overlap matrix
  for (int i = 0; i < no; i++)
  for (int j = 0; j < no; j++) {
    int p, q;
    S val;
    ifs >> p >> q >> val;
    ovi[p*no+q] = val;
  }

  // read the one electron integrals
  for (int i = 0; i < no; i++)
  for (int j = 0; j < no; j++) {
    int p, q;
    S val;
    ifs >> p >> q >> val;
    oei[p*no+q] = val;
  }

  // read the two electron integrals
  {
    std::fill(tei, tei + formic::tei_compact_array_length(no), formic::zero(S()));
    int p, q, r, s;
    S val;
    while ( ifs >> p >> q >> r >> s >> val )
      tei[formic::tei_compound_index(p,q,r,s)] = val;
  }

  // close the input file
  ifs.close();

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns a pointer to the dense two electron integrals in 1122 storage
/// 
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> const S * formic::Hamiltonian<S>::dense_tei_ptr() {

  // if the integrals are not present in dense storage, prepare them
  if ( _dense_tei.size() != _norb*_norb*_norb*_norb ) {
    _dense_tei.reset(_norb*_norb*_norb*_norb);
    formic::tei_compact_to_dense(_norb, &_tei.at(0), &_dense_tei.at(0));
    //formic::tei_dense_to_compact(_norb, &_dense_tei.at(0), &_tei.at(0));
    //formic::tei_compact_to_dense(_norb, &_tei.at(0), &_dense_tei.at(0));
  }

  return &_dense_tei.at(0);

}

template const double * formic::Hamiltonian<double>::dense_tei_ptr();
template const std::complex<double> * formic::Hamiltonian<std::complex<double> >::dense_tei_ptr();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   creates text files with the one and two electron integrals in dmrg input format
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class S> void formic::Hamiltonian<S>::write_dmrg_ints(const formic::InputBase & userinp) const {

  // get the file names for the one and two electron integrals
  const std::string int1_file = "dmrg_int1.txt";
  const std::string int2_file = "dmrg_int2.txt";

  // open the one electron integral file
  std::ofstream ofs_int1(int1_file.c_str());

  // check that the file is open
  if (!ofs_int1.is_open())
    throw formic::Exception("failed to open file \"%s\" in formic::Hamiltonian::write_dmrg_ints") % int1_file;

  // write the number of orbitals
  ofs_int1 << _norb << std::endl;

  // write the one electron integrals
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q < _norb; q++)
    ofs_int1 << boost::format("%4i  %4i  %60s\n") % p % q % formic::format_number("%22.12e", this->get_oei(p, q));

  // close the one electron integral file
  ofs_int1.close();

  // open the two electron integral file
  std::ofstream ofs_int2(int2_file.c_str());

  // check that the file is open
  if (!ofs_int2.is_open())
    throw formic::Exception("failed to open file \"%s\" in formic::Hamiltonian::write_dmrg_ints") % int2_file;

  // write the number of orbitals
  ofs_int2 << _norb << std::endl;

  // write the two electron integrals
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q < _norb; q++)
  for (int r = 0; r < _norb; r++)
  for (int s = 0; s < _norb; s++)
    ofs_int2 << boost::format("%4i  %4i  %4i  %4i  %60s\n") % p % r % q % s % formic::format_number("%22.12e", this->get_tei(p, q, r, s));

  // close the two electron integral file
  ofs_int2.close();

}

// explicitly instantiate the function
template void formic::Hamiltonian< double               >::write_dmrg_ints(const formic::InputBase &) const;
template void formic::Hamiltonian< std::complex<double> >::write_dmrg_ints(const formic::InputBase &) const;

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   writes the integrals to a file in the text format used by formic
///
/// \param[in]   userinp   user input options
/// \param[in]   filename  the name of the file to write to
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::Hamiltonian<S>::write_text_ints(const formic::InputBase & userinp, const std::string & filename) const {

  // open the file
  std::ofstream ofs(filename.c_str());

  // check that the file is open
  if (!ofs.is_open())
    throw formic::Exception("failed to open file \"%s\" in formic::Hamiltonian::write_text_ints") % filename;

  // write the number of orbitals
  ofs << _norb << std::endl << std::endl;

//  // write the number of irreps
//  ofs << userinp.get<std::vector<int> >("opirr").size() << std::endl << std::endl;
//
//  // write the orbital irreps
//  for (int i = 0; i < userinp.get<std::vector<int> >("orb_irrep").size(); i++)
//    ofs << userinp.get<std::vector<int> >("orb_irrep").at(i) << " ";
//  ofs << std::endl << std::endl;

  // write nuclear repulsion energy
  ofs << boost::format("%20.12f") % _e_cor_nuc << std::endl << std::endl;

  // write overlap integrals
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++)
    ofs << boost::format(" %4i %4i %20.12f") % i % j % ( i == j ? 1.0 : 0.0 ) << std::endl;
  ofs << std::endl;

  // write one electron integrals
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++)
    ofs << boost::format(" %4i %4i %20.12f") % i % j % _oei[i*_norb+j] << std::endl;
  ofs << std::endl;

  // write two electron integrals
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j <= i;    j++)
  for (int k = 0; k < _norb; k++)
  for (int l = 0; l <= k;    l++)
  if ( this->get_tei(i,j,k,l) != formic::zero(S()) )
    ofs << boost::format(" %4i %4i %4i %4i %20.12f") % i % j % k % l % this->get_tei(i,j,k,l) << std::endl;
  ofs << std::endl;

  // close the file
  ofs.close();

}

template void formic::Hamiltonian<double>::write_text_ints(const formic::InputBase &, const std::string &) const;
template void formic::Hamiltonian<std::complex<double> >::write_text_ints(const formic::InputBase &, const std::string &) const;

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   prepares the Hubbard model Hamiltonian for fermions or spinless fermions
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class S> void formic::Hamiltonian<S>::initialize_hubbard(const formic::InputBase & userinp) {

  throw formic::Exception("not yet implemented");

}

template void formic::Hamiltonian< double               >::initialize_hubbard(const formic::InputBase &);
template void formic::Hamiltonian< std::complex<double> >::initialize_hubbard(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the modified one electron integrals
///
///   The quantum chemistry electronic hamiltonian can be written as
///
///   \f[ H = \sum_{i j} h_{i j} a^\dag_i a_j + \frac{1}{2} \sum_{ p q r s} [pq|rs] a^\dag_p a^\dag_r a_s a_q \f]
///
///   where \f$h_{i j}\f$ are the one electron integrals stored in _oei and
///   \f$[pq|rs]\f$ are the two electron integrals stored in _tei.
///
///   The local energy function in quantum monte carlo is more easily computed if we
///   rearrange the Hamiltonian as 
///
///   \f[ H = \sum_{i j} f_{i j} a^\dag_i a_j + \frac{1}{2} \sum_{ p q r s} [pq|rs] a^\dag_p a_q a^\dag_r a_s \f]
///
///   where we have commuted \f$a_q\f$ to be next to \f$a^\dag_p\f$ and absorbed the
///   result of the delta function with \f$a^\dag_r\f$ into the one electron integrals,
///   creating the modified one electron integrals
///
///   \f[  f_{i j} = h_{i j} - \frac{1}{2} \sum_p  [ip|pj]  \f]
///
///   These modified one electron integrals are stored in _mei.
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class S> void formic::Hamiltonian<S>::initialize_mei() {

  // modify the one electron integrals
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++) {
    _mei[i*_norb+j] = _oei[i*_norb+j];
    for (int p = 0; p < _norb; p++)
      _mei[i*_norb+j] -= 0.5 * this->get_tei(i,p,p,j);
  }

}

template void formic::Hamiltonian< double               >::initialize_mei();
template void formic::Hamiltonian< std::complex<double> >::initialize_mei();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the cholesky decomposition of the two electron integrals
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::Hamiltonian<S>::compute_cholesky() {

    formic::start_timer("cholesky decomposition");
    formic::cholesky_decomp(*this, *this);
    formic::stop_timer("cholesky decomposition");
    formic::of << boost::format("cholesky threshold of %.12f resulted in %i cholesky vectors")
                   % this->chol_thresh()
                   % _n_cholesky
                << std::endl
                << std::endl;

}

template void formic::Hamiltonian<double>::compute_cholesky();
template void formic::Hamiltonian<std::complex<double> >::compute_cholesky();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   loads the quantum chemistry integrals and computes their cholesky decomposition
///
/// \param[in]   userinp  user input options
/// \param[in]   source   a string to specify where the integrals should be read from
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::Hamiltonian<S>::initialize_integrals(const formic::InputBase & userinp,
                                                                     const std::string & source,
                                                                     const std::string & filename) {

  // print some info
  if ( formic::mpi::rank() == 0 ) {
    formic::of << boost::format("   ==================================================") << std::endl;
    formic::of << boost::format("   ======   Initializing Integrals From %4s   ======") % source << std::endl;
    formic::of << boost::format("   ==================================================") << std::endl;
    formic::of << std::endl;
  }

  // set the number of orbitals
  _norb = userinp.get<int>("nao");

  // set the cholesky decomposition threshold
  _chol_thresh = userinp.get<double>("chol_thresh");

  // reset the cholesky vectors
  _n_cholesky = 0;
  _cholesky_vecs.clear();
  _cholesky_vecs_t.clear();

  // empty the dense two electron integrals (they will be rebuilt the next time they are needed)
  _dense_tei.clear();

  // allocate storage for the integrals
  _ovi.assign(_norb*_norb, formic::zero(S()));
  _oei.assign(_norb*_norb, formic::zero(S()));
  _mei.assign(_norb*_norb, formic::zero(S()));
  _tei.assign(formic::tei_compact_array_length(_norb), formic::zero(S()));

  // set the initial orbital coefficient matrix to the identity
  _orb_coeff = formic::identity_matrix<S,int>(_norb);

  // read in integrals on root node
  if ( formic::mpi::rank() == 0 ) {

    // read integrals on root process
    if (source == "psi3")
      formic::get_psi3_ints(userinp, _e_cor_nuc, &_ovi.at(0), &_oei.at(0), &_tei.at(0));
    //else if (source == "molpro")
    //  formic::get_molpro_ints(userinp, _e_cor_nuc, &_oei.at(0), &_tei.at(0));
    else if (source == "text")
      formic::get_text_ints(userinp, filename, _e_cor_nuc, &_ovi.at(0), &_oei.at(0), &_tei.at(0));
    else
      throw formic::Exception("unknown integral source \"%s\" in formic::Hamiltonian::initialize_integrals") % source;

    // initialize the modified one electron integrals
    this->initialize_mei();

    // compute the cholesky decomposition of the two electron integrals
    this->compute_cholesky();

    // write initial hamiltonian to disk
    this->to_archive("fqmc_init_hamiltonian.bin");

//    // compute orthonormal orbitals
//    formic::orthogonalize_orbs_by_overlap(_norb, _ovi, _oei, _tei);

  }

  // read the hamiltonian from disk (this also broadcasts it to all processes)
  this->from_archive("fqmc_init_hamiltonian.bin");

//  // transform orbitals to orthonormal S^(-1/2) basis and send them to all processes
//  std::vector<S> s_neg_half(_ovi);
//  if ( formic::mpi::rank() == 0 )
//    formic::matrix_neg_half(_norb, &s_neg_half.at(0));
//  this->transform_orbitals(&s_neg_half.at(0));

  // print initial integrals
  if ( formic::mpi::rank() == 0 ) {
    this->print_oei(_norb, _oei, "initial one electron integrals:");
    this->print_oei(_norb, _ovi, "initial overlap integrals:");
    this->print_tei_diag(_norb, _tei, "[pp|qq] elements of initial two electron integrals:");
  }

//  // print nuclear repulsion energy
//  formic::of << boost::format("enuc = %20.12f") % _e_cor_nuc << std::endl;
//  formic::of << std::endl;
//
//  // print one electron integrals
//  formic::of << boost::format("one electron integrals:") << std::endl;
//  formic::of << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//    formic::of << boost::format(" %4i %4i %20.12f") % i % j % _oei[i*_norb+j] << std::endl;
//  formic::of << std::endl;
//
//  // print two electron integrals
//  formic::of << boost::format("two electron integrals:") << std::endl;
//  formic::of << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//  for (int k = 0; k < _norb; k++)
//  for (int l = 0; l < _norb; l++)
//    formic::of << boost::format(" %4i %4i %4i %4i %20.12f") % i % j % k % l % this->get_tei(i,j,k,l) << std::endl;
//  formic::of << std::endl;
//
//  // print two electron integral diagonal
//  std::cout << boost::format("diagonal two electron integrals:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//    std::cout << boost::format(" %4i %4i %20.12f") % i % j % this->get_tei(i,j,i,j) << std::endl;
//  std::cout << std::endl;
//
//  for (int i = 0; i < _cholesky_vecs.size(); i++) {
//    if (i % (_norb*_norb) == 0)
//      std::cout << std::endl;
//    std::cout << boost::format("cholesky %6i %20.12f") % i % _cholesky_vecs.at(i) << std::endl;
//  }
//
//  // print two electron integral error
//  double error = 0.0;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//  for (int k = 0; k < _norb; k++)
//  for (int l = 0; l < _norb; l++) {
//    S x = this->get_tei(i,j,k,l);
//    for (int p = 0; p < _n_cholesky; p++)
//      x -= _cholesky_vecs.at(_norb*_norb*p + _norb*i + j) * _cholesky_vecs.at(_norb*_norb*p + _norb*k + l);
//    if (std::abs(x) > this->chol_thresh())
//      std::cout << boost::format("error in integral %4i %4i %4i %4i %20.12f") % i % j % k % l % this->get_tei(i,j,k,l)
//                << std::endl;
//    error += formic::square_norm(x);
//  }
//  std::cout << std::endl;
//  error = std::sqrt(error);
//  std::cout << boost::format("total error in two electron integrals is %22.12e") % error << std::endl;

}

template void formic::Hamiltonian< double               >::initialize_integrals(const formic::InputBase &, const std::string &, const std::string &);
template void formic::Hamiltonian< std::complex<double> >::initialize_integrals(const formic::InputBase &, const std::string &, const std::string &);


template<class T> void formic::Hamiltonian<T>::to_archive(formic::Archive & ar) const {

  ar << this->_norb;
  ar << this->_e_cor_nuc;
  ar << this->_orb_coeff;
  ar << this->_ovi;
  ar << this->_oei;
  ar << this->_mei;
  ar << this->_tei;
  ar << this->_n_cholesky;
  ar << this->_cholesky_vecs;
  //ar << this->_cholesky_vecs_t;

}

template<class T> void formic::Hamiltonian<T>::to_archive(const std::string & filename) const {

  // only the root process writes to disk
  if (formic::mpi::rank() == 0) {
    formic::BinaryFileArchive ar(filename, "wt");
    ar << *this;
  }

}

template<class T> void formic::Hamiltonian<T>::from_archive(formic::Archive & ar) {

  // empty the dense two electron integrals (they will be rebuilt the next time they are needed)
  _dense_tei.clear();

  // read hamiltonian data from the archive
  formic::mpi::read_and_bcast(ar, this->_norb);
  formic::mpi::read_and_bcast(ar, this->_e_cor_nuc);
  formic::mpi::read_and_bcast(ar, this->_orb_coeff);
  formic::mpi::read_and_bcast(ar, this->_ovi);
  formic::mpi::read_and_bcast(ar, this->_oei);
  formic::mpi::read_and_bcast(ar, this->_mei);
  formic::mpi::read_and_bcast(ar, this->_tei);
  formic::mpi::read_and_bcast(ar, this->_n_cholesky);
  formic::mpi::read_and_bcast(ar, this->_cholesky_vecs);
  //formic::mpi::read_and_bcast(ar, this->_cholesky_vecs_t);

  // compute transpose of cholesky vectors
  _cholesky_vecs_t.resize(_cholesky_vecs.size());
  for (int j = 0; j < _n_cholesky; j++)
    formic::xcopy(_norb*_norb, &_cholesky_vecs.at(_norb*_norb*j), 1, &_cholesky_vecs_t.at(j), _n_cholesky);

}

template<class T> void formic::Hamiltonian<T>::from_archive(const std::string & filename) {

  // root process reads from disk
  if (formic::mpi::rank() == 0) {

    formic::BinaryFileArchive ar(filename, "r");
    ar >> *this;

  // other processes get their data from root and so we just pass them an empty archive
  } else {

    formic::BinaryArchive ar;
    ar >> *this;

  }

}

template<class T, class S> formic::BinaryArchiveTemplate<T> & operator<<(formic::BinaryArchiveTemplate<T> & ar, const formic::Hamiltonian<S> & ham) {

  // root process writes to disk
  if (formic::mpi::rank() == 0)
    ham.to_archive(ar);

  // return the archive
  return ar;

}

template formic::BinaryArchiveTemplate<std::stringstream> & operator<<(formic::BinaryArchiveTemplate<std::stringstream> &,
                                                                       const formic::Hamiltonian<double> &);
template formic::BinaryArchiveTemplate<std::stringstream> & operator<<(formic::BinaryArchiveTemplate<std::stringstream> &,
                                                                       const formic::Hamiltonian<std::complex<double> > &);
template formic::BinaryArchiveTemplate<std::fstream> & operator<<(formic::BinaryArchiveTemplate<std::fstream> &,
                                                                  const formic::Hamiltonian<double> &);
template formic::BinaryArchiveTemplate<std::fstream> & operator<<(formic::BinaryArchiveTemplate<std::fstream> &,
                                                                  const formic::Hamiltonian<std::complex<double> > &);

template<class T, class S> formic::BinaryArchiveTemplate<T> & operator>>(formic::BinaryArchiveTemplate<T> & ar, formic::Hamiltonian<S> & ham) {

  // root process reads from disk
  if (formic::mpi::rank() == 0) {

    ham.from_archive(ar);

  // other processes get their data from root and so we just pass them an empty archive
  } else {

    formic::BinaryArchive empty;
    ham.from_archive(empty);

  }

  // return the archive
  return ar;

}

template formic::BinaryArchiveTemplate<std::stringstream> & operator>>(formic::BinaryArchiveTemplate<std::stringstream> &,
                                                                       formic::Hamiltonian<double> &);
template formic::BinaryArchiveTemplate<std::stringstream> & operator>>(formic::BinaryArchiveTemplate<std::stringstream> &,
                                                                       formic::Hamiltonian<std::complex<double> > &);
template formic::BinaryArchiveTemplate<std::fstream> & operator>>(formic::BinaryArchiveTemplate<std::fstream> &,
                                                                  formic::Hamiltonian<double> &);
template formic::BinaryArchiveTemplate<std::fstream> & operator>>(formic::BinaryArchiveTemplate<std::fstream> &,
                                                                  formic::Hamiltonian<std::complex<double> > &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Get the [pq,pq] two electron integrals.
///
/// \param[out]     pqpq_ints   on exit, the [pq,pq] integrals
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class T> void formic::Hamiltonian<T>::get_pqpq_ints(std::vector<T> & pqpq_ints) const {

  pqpq_ints.resize(_norb*_norb);
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q < _norb; q++) {
    pqpq_ints[_norb*p+q] = formic::xdotc(_n_cholesky, &_cholesky_vecs_t[_n_cholesky * ( _norb * p + q )], 1,
                                                      &_cholesky_vecs_t[_n_cholesky * ( _norb * p + q )], 1);
    //pqpq_ints[_norb*p+q] = formic::zero<T>();
    //for (int u = 0; u < _n_cholesky; u++) {
    //  const T x = _cholesky_vecs_t[ _n_cholesky * ( _norb * p + q ) + u ];
    //  pqpq_ints[_norb*p+q] += x * x;
    //}
  }

}

template void formic::Hamiltonian<double>::get_pqpq_ints(std::vector<double> &) const;
template void formic::Hamiltonian<std::complex<double> >::get_pqpq_ints(std::vector<std::complex<double> > &) const;

template<class T> void formic::Hamiltonian<T>::print_oei(const int n, const std::vector<T> & mat, const std::string & title) {

  formic::of << title << std::endl;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++)
      formic::of << boost::format(" %10.6f") % mat.at(i*n+j);
    formic::of << std::endl;
  }
  formic::of << std::endl;

}

template<class T> void formic::Hamiltonian<T>::print_tei_diag(const int n, const std::vector<T> & mat, const std::string & title) {

  formic::of << title << std::endl;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++)
      formic::of << boost::format(" %10.6f") % mat.at(formic::tei_compound_index(i,i,j,j));
    formic::of << std::endl;
  }
  formic::of << std::endl;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   converts from the AO 1- and 2-electron integrals to the S^(-1/2) orthogonalized basis
///          1- and 2- electron integrals
///
/// \param[in]      n        the number of ao orbitals
/// \param[in]      overlap  the ao overlap matrix, size n*n
/// \param[in,out]  int1     on entry, the ao one electron integrals, size n*n
///                          on exit, the one electron integrals in the new basis
/// \param[in,out]  int2     on entry, the ao two electron integrals,
///                          size  n*(n+1)/2 * ( n*(n+1)/2 + 1 ) / 2
///                          on exit, the two electron integrals in the new basis
///
///////////////////////////////////////////////////////////////////////////////////////////////////

void formic::orthogonalize_orbs_by_overlap(const int n,
                                           const std::vector<double> & overlap,
                                           std::vector<double> & int1,
                                           std::vector<double> & int2)
{

  // allocate some work space
  std::vector<double> intermed_0(n*n, 0.0);
  std::vector<double> intermed_1(n*n, 0.0);
  std::vector<double> big2_0(n*n*n*n, 0.0);
  std::vector<double> big2_1(n*n*n*n, 0.0);

//  // preparation for diagonalization of overlap matrix
//  int info;
//  const int lwork = 5*n;
//  std::vector<double> work(lwork, 0.0);
//  std::vector<double> rwork(3*n, 0.0);
//  std::vector<double> eval(n, 0.0);
//
//  // diagonalize the overlap matrix, after which the eigenvectors will be in the rows the matrix
//  std::vector<double> s_neg_half(overlap);
//  formic::xsyev('V', 'U', n, &s_neg_half.at(0), n, &eval.at(0), &work.at(0), lwork, &rwork.at(0), info);
//  if (info != 0)
//    throw formic::Exception("formic::xsyev failed with error code %i in formic:::orthogonalize_orbs_by_overlap") % info;
//
//  // compute the S^(-1/2) matrix
//  for (int i = 0; i < n; i++) {
//    const double x = 1.0 / std::sqrt(eval[i]);
//    for (int j = 0; j < n; j++) {
//      intermed_0[i*n+j] = s_neg_half[j*n+i];
//      intermed_1[i*n+j] = s_neg_half[i*n+j] * x;
//    }
//  }
//  formic::xgemm('N', 'N', n, n, n, 1.0, &intermed_1.at(0), n, &intermed_0.at(0), n, 0.0, &s_neg_half.at(0), n);

  // compute the S^(-1/2) matrix
  std::vector<double> s_neg_half(overlap);
  formic::matrix_neg_half(n, &s_neg_half.at(0));

  // print the orthonormal orbital coefficients
  formic::Hamiltonian<double>::print_oei(n, s_neg_half, "orthonormal basis orbitals coefficients for S^(-1/2) basis:");

  // convert the one electron integrals
  formic::xgemm('N', 'N', n, n, n, 1.0, &s_neg_half.at(0), n, &int1.at(0), n, 0.0, &intermed_0.at(0), n);
  formic::xgemm('N', 'N', n, n, n, 1.0, &intermed_0.at(0), n, &s_neg_half.at(0), n, 0.0, &int1.at(0), n);

  // print the one electron integrals
  formic::Hamiltonian<double>::print_oei(n, int1, "one electron integrals for S^(-1/2) basis:");

  // convert the two electron integrals

  for (int p = 0; p < n; p++)
  for (int q = 0; q < n; q++)
  for (int r = 0; r < n; r++)
  for (int s = 0; s < n; s++)
    big2_0[n*n*n*p+n*n*q+n*r+s] = int2[formic::tei_compound_index(p,q,r,s)];

  struct ConvertAndShift {
    static void foo(const int n, const std::vector<double> & s_neg_half, double * p0, double * p1) {
      formic::xgemm('N', 'N', n*n*n, n, n, 1.0, p0, n*n*n, &s_neg_half.at(0), n, 0.0, p1, n*n*n);
      for (int p = 0; p < n; p++)
      for (int q = 0; q < n; q++)
      for (int r = 0; r < n; r++)
      for (int s = 0; s < n; s++)
        p0[n*n*n*p+n*n*q+n*r+s] = p1[n*n*n*q+n*n*r+n*s+p];
    }
  };
  ConvertAndShift::foo(n, s_neg_half, &big2_0.at(0), &big2_1.at(0));
  ConvertAndShift::foo(n, s_neg_half, &big2_0.at(0), &big2_1.at(0));
  ConvertAndShift::foo(n, s_neg_half, &big2_0.at(0), &big2_1.at(0));
  ConvertAndShift::foo(n, s_neg_half, &big2_0.at(0), &big2_1.at(0));

  for (int p = 0; p < n; p++)
  for (int q = 0; q < n; q++)
  for (int r = 0; r < n; r++)
  for (int s = 0; s < n; s++)
    int2[formic::tei_compound_index(p,q,r,s)] = big2_0[n*n*n*p+n*n*q+n*r+s];

  // print the [pp|qq] two electron integrals
  formic::Hamiltonian<double>::print_tei_diag(n, int2, "[pp|qq] two electron integrals for S^(-1/2) basis:");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   converts from the AO 1- and 2-electron integrals to the S^(-1/2) orthogonalized basis
///          1- and 2- electron integrals
///
/// \param[in]      n        the number of ao orbitals
/// \param[in]      overlap  the ao overlap matrix, size n*n
/// \param[in,out]  int1     on entry, the ao one electron integrals, size n*n
///                          on exit, the one electron integrals in the new basis
/// \param[in,out]  int2     on entry, the ao two electron integrals,
///                          size  n*(n+1)/2 * ( n*(n+1)/2 + 1 ) / 2
///                          on exit, the two electron integrals in the new basis
///
///////////////////////////////////////////////////////////////////////////////////////////////////

void formic::orthogonalize_orbs_by_overlap(const int n,
                                           const std::vector<std::complex<double> > & overlap,
                                           std::vector<std::complex<double> > & int1,
                                           std::vector<std::complex<double> > & int2)
{

  throw formic::Exception("overlap orthogonalization not yet implemented for complex numbers");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Reorders the atomic orbitals
///
/// \param[in]      userinp     the user's input options
/// \param[in]      order       the new order for the orbitals, size _norb
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class T> void formic::Hamiltonian<T>::reorder_orbitals(const formic::InputBase & userinp,
                                                                const std::vector<int> & order) {

  //throw formic::Exception("formic::Hamiltonian::reorder_orbitals has not been updated to work with the new orbital transformation system");

  // make sure new order is sane
  if ( formic::mpi::rank() == 0 ) {
    if ( order.size() != _norb )
      throw formic::Exception("new order vector in formic::Hamiltonian::reorder_orbitals must have length equal to the number of orbitals");
    for (int p = 0; p < _norb; p++) {
      if ( order[p] < 0 || order[p] >= _norb )
        throw formic::Exception("bad orbital index \"%i\" in new order vector in formic::Hamiltonian::reorder_orbitals") % order[p];
      for (int q = 0; q < p; q++)
        if ( order[p] == order[q] )
          throw formic::Exception("repeated orbital index \"%i\" in new order vector in formic::Hamiltonian::reorder_orbitals") % order[p];
    }
  }

  // prepare orbital transformation matrix for the reordering
  std::vector<T> transform_matrix(_norb*_norb, formic::zero(T()));
  if ( formic::mpi::rank() == 0 )
    for (int i = 0; i < _norb; i++)
      transform_matrix.at(i+order.at(i)*_norb) = formic::unity(T());

  this->transform_orbitals(&transform_matrix.at(0));

//  // root node does the reordering
//  if ( formic::mpi::rank() == 0 ) {
//
//    // make sure new order is sane
//    if ( order.size() != _norb )
//      throw formic::Exception("new order vector in formic::Hamiltonian::reorder_orbitals must have length equal to the number of orbitals");
//    for (int p = 0; p < _norb; p++) {
//      if ( order[p] < 0 || order[p] >= _norb )
//        throw formic::Exception("bad orbital index \"%i\" in new order vector in formic::Hamiltonian::reorder_orbitals") % order[p];
//      for (int q = 0; q < p; q++)
//        if ( order[p] == order[q] )
//          throw formic::Exception("repeated orbital index \"%i\" in new order vector in formic::Hamiltonian::reorder_orbitals") % order[p];
//    }
//
//    const int noei = _norb*_norb;
//    const int ntei = formic::tei_compact_array_length(_norb);
//
//    // get temporary workspace
//    std::vector<T> work(std::max(noei, ntei), formic::zero(T()));
//
//    // reorder overlap integrals
//    if ( _ovi.size() == noei ) {
//      std::vector<T> & ref = _ovi;
//      for (int p = 0; p < _norb; p++)
//        for (int q = 0; q < _norb; q++)
//          work[ order[p] * _norb + order[q] ] = ref[ p * _norb+ q ];
//      std::copy(work.begin(), work.begin() + ref.size(), ref.begin());
//    }
//
//    // reorder one electron integrals
//    if ( _oei.size() == noei ) {
//      std::vector<T> & ref = _oei;
//      for (int p = 0; p < _norb; p++)
//        for (int q = 0; q < _norb; q++)
//          work[ order[p] * _norb + order[q] ] = ref[ p * _norb+ q ];
//      std::copy(work.begin(), work.begin() + ref.size(), ref.begin());
//    }
//
//    // reorder modified one electron integrals
//    if ( _mei.size() == noei ) {
//      std::vector<T> & ref = _mei;
//      for (int p = 0; p < _norb; p++)
//        for (int q = 0; q < _norb; q++)
//          work[ order[p] * _norb + order[q] ] = ref[ p * _norb+ q ];
//      std::copy(work.begin(), work.begin() + ref.size(), ref.begin());
//    }
//
//    // reorder two electron integrals
//    if ( _tei.size() == ntei ) {
//      for (int p = 0; p < _norb; p++)
//      for (int q = 0; q < _norb; q++)
//      for (int r = 0; r < _norb; r++)
//      for (int s = 0; s < _norb; s++)
//        work[ formic::tei_compound_index( order[p], order[q], order[r], order[s] ) ] = _tei[formic::tei_compound_index(p,q,r,s)];
//      std::copy(work.begin(), work.begin() + _tei.size(), _tei.begin());
//      // re-compute cholesky decomposition
//      this->compute_cholesky();
//    }
//
//    // print reordered integrals
//    this->print_oei(_norb, _oei, "reordered one electron integrals:");
//    this->print_tei_diag(_norb, _tei, "[pp|qq] elements of reordered two electron integrals:");
//
//    // write hamiltonian to disk
//    this->to_archive("hamiltonian_dump.bin");
//
//  }
//
//  // read the hamiltonian from disk (this also broadcasts it to all processes)
//  this->from_archive("hamiltonian_dump.bin");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "Psi3Ints" command, which loads integrals from Psi3
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////

static void formic_Psi3Ints_command_function(formic::Archive & arch, formic::AllData & data) {

  // load orbs per irrep and orbital irreps
  {

    // on the root process, load the number of orbitals in each irrep
    std::vector<int> opirr;
    if ( formic::mpi::rank() == 0 )
      formic::get_psi3_orbs_per_irrep(data.userinp, opirr);

    // run the Option command to update the input options with the number of orbs per irrep
    formic::run_command("Option", formic::TextArchive() << "opirr" << opirr, data);

    // run the Option command to update the input options with the orbital irreps
    const int nao = data.userinp.get<int>("nao");
    std::vector<int> orb_irrep(nao, 0);
    if ( formic::mpi::rank() == 0 )
      for (int i = 0, p = 0; i < opirr.size(); i++)
        for (int j = 0; j < opirr.at(i); j++)
          orb_irrep.at(p++) = i;
    formic::run_command("Option", formic::TextArchive() << "orb_irrep" << orb_irrep, data);

  }

  // print orbs per irrep
  if ( formic::mpi::rank() == 0 ) {
    formic::of << boost::format("Psi3 orbs per irrep:  ");
    for (int i = 0; i < data.userinp.get<std::vector<int> >("opirr").size(); i++)
      formic::of << boost::format(" %4i") % data.userinp.get<std::vector<int> >("opirr").at(i);
    formic::of << std::endl;
    formic::of << std::endl;
  }

  // print orbital irreps
  if ( formic::mpi::rank() == 0 ) {
    formic::of << boost::format("Psi3 orbital irreps:") << std::endl;
    for (int i = 0; i < data.userinp.get<std::vector<int> >("orb_irrep").size(); i++)
      formic::of << boost::format(" %4i  %i") % i % data.userinp.get<std::vector<int> >("orb_irrep").at(i) << std::endl;
    formic::of << std::endl;
  }

  // load the integrals from psi3
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    data.hd.initialize_integrals(data.userinp, "psi3", "");
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    data.hc.initialize_integrals(data.userinp, "psi3", "");
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_Psi3Ints_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "TextInts" command, which loads integrals from a text file
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////

static void formic_TextInts_command_function(formic::Archive & arch, formic::AllData & data) {

  // get the filename
  std::string filename;
  formic::mpi::read_and_bcast(arch, filename, "failed to read filename for command \"TextInts\"");

//  // read in orbitals per irrep on root process
//  std::vector<int> opirr;
//  std::vector<int> orb_irrep;
//  if (formic::mpi::rank() == 0) {
//
//    // open the file
//    std::ifstream ifs(filename.c_str());
//    if (!ifs.is_open())
//      throw formic::Exception("failed to open file \"%s\" in formic_TextInts_command_function") % filename;
//
//    // read the number of orbitals
//    int norb = 0;
//    ifs >> norb;
//    if ( norb != data.userinp.get<int>("nao") )
//      throw formic::Exception("wrong number of orbitals \"%i\" in formic_TextInts_command_function") % norb;
//
//    // read the number of irreps
//    int nirr = 0;
//    ifs >> nirr;
//
//    // read the orbital of each irrep
//    opirr.assign(nirr, 0);
//    orb_irrep.assign(norb, 0);
//    for (int i = 0; i < norb; i++) {
//      ifs >> orb_irrep.at(i);
//      opirr.at(orb_irrep.at(i)) += 1;
//    }
//
//    // close the file
//    ifs.close();
//
//  }
//
//  // run the Option command to update the input options with the number of orbs per irrep
//  formic::run_command("Option", formic::TextArchive() << "opirr" << opirr, data);
//
//  // run the Option command to update the input options with the orbital irreps
//  formic::run_command("Option", formic::TextArchive() << "orb_irrep" << orb_irrep, data);

  // load the integrals from the file
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    data.hd.initialize_integrals(data.userinp, "text", filename);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    data.hc.initialize_integrals(data.userinp, "text", filename);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_TextInts_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "ReorderOrbitals" command
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_ReorderOrbitals_command_function(formic::Archive & arch, formic::AllData & data) {

  std::vector<int> order;
  formic::mpi::read_and_bcast(arch, order, "failed to read order in formic_ReorderOrbitals_command_function");

  // reorder orbital irreps
  {

    // run the Option command to update the input options with the orbital irreps
    const int nao = data.userinp.get<int>("nao");
    std::vector<int> orb_irrep(nao, 0);
    if ( formic::mpi::rank() == 0 )
      for (int i = 0; i < nao; i++)
        orb_irrep.at(i) = data.userinp.get<std::vector<int> >("orb_irrep").at(order.at(i));
    formic::run_command("Option", formic::TextArchive() << "orb_irrep" << orb_irrep, data);

  }

  // print orbital irreps
  if ( formic::mpi::rank() == 0 ) {
    formic::of << boost::format("reordered orbital irreps:") << std::endl;
    for (int i = 0; i < data.userinp.get<std::vector<int> >("orb_irrep").size(); i++)
      formic::of << boost::format(" %4i  %i") % i % data.userinp.get<std::vector<int> >("orb_irrep").at(i) << std::endl;
    formic::of << std::endl;
  }

  if (data.userinp.get<std::string>("scalar_type") == "real") {
    data.hd.reorder_orbitals(data.userinp, order);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    data.hc.reorder_orbitals(data.userinp, order);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_ReorderOrbitals_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "WriteTextInts" command, which write's the integrals stored
///          in the Hamiltonian to a text file that can be read in with the TextInts command.
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_WriteTextInts_command_function(formic::Archive & arch, formic::AllData & data) {

  // get the filename to write the integrals to
  std::string filename;
  formic::mpi::read_and_bcast(arch, filename, "failed to read filename in formic_WriteTextInts_command_function");

  // write out the integrals on the root process
  if ( formic::mpi::rank() == 0 ) {
    if (data.userinp.get<std::string>("scalar_type") == "real") {
      data.hd.write_text_ints(data.userinp, filename);
    } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
      data.hc.write_text_ints(data.userinp, filename);
    } else 
      throw formic::Exception("Unknown scalar type \"%s\" in formic_WriteTextInts_command_function.") % data.userinp.get<std::string>("scalar_type");
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "ReadBinHam" command, which reads the hamiltonian from
///          a binary archive file.
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_ReadBinHam_command_function(formic::Archive & arch, formic::AllData & data) {

  // get the name of the file to read from
  std::string filename;
  formic::mpi::read_and_bcast(arch, filename, "failed to read filename in formic_ReadBinHam_command_function");

  // say what we are doing
  if ( formic::mpi::rank() == 0 )
    formic::of << boost::format("reading hamiltonian from binary file \"%s\"") % filename << std::endl << std::endl;

  // read in the hamiltonian
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    data.hd.from_archive(filename);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    data.hc.from_archive(filename);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_ReadBinHam_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function used by user commands that set or rotate the orbital basis
///
/// \param[in]     userinp      user input options
/// \param[in]     ham          the hamiltonian
/// \param[in]     wfn          wave function container
/// \param[in]     repl_inds    indices of orbitals to replace with new/rotated orbitals
/// \param[in]     new_inds     indices of orbitals to use for the replacement/rotation
/// \param[in]     oc_name      name of the data object to get new coefficients from
/// \param[in]     rotate       whether we are rotating (if true) or setting (if false)
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
static void formic_set_or_rotate_basis_helper_function(const formic::InputBase & userinp,
                                                       formic::Hamiltonian<S> & ham,
                                                       formic::Wavefunction<S> & wfn,
                                                       const std::vector<int> & repl_inds,
                                                       const std::vector<int> & new_inds,
                                                       const std::string & oc_name,
                                                       const bool rotate) {

  // get number of orbitals
  const int norb = userinp.get<int>("nao");

  // check that new and replacement index vectors are the same length
  if ( repl_inds.size() != new_inds.size() )
    throw formic::Exception("new and replacement index vector lengths do not match in formic_set_or_rotate_basis_helper_function");

  // get the coefficients from which we will build the replacement/rotation
  std::vector<S> & oc_vec = wfn.get_wf(oc_name)->data_vec();

  // check that the coefficient data is the correct length
  if ( oc_vec.size() != norb*norb )
    throw formic::Exception("coefficient data vector length is incorrect in formic_set_or_rotate_basis_helper_function");

  // get initial replacement or rotation matrix
  formic::Matrix<S,int> rot_repl_mat = ( rotate ? formic::identity_matrix<S,int>(norb) : ham.orb_coeff().clone() );

  // update the replacement/rotation matrix with the new coefficients
  for (int i = 0; i < repl_inds.size(); i++)
    for (int j = 0; j < norb; j++)
      rot_repl_mat.at(repl_inds.at(i), j) = oc_vec.at(j + norb * new_inds.at(i));

  // execute the basis set transformation
  ham.transform_orbitals(&rot_repl_mat.at(0,0), rotate);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "RotateBasisOrbsSelect" command, which replaces a select set
///          of the current basis orbitals with chosen rotations of the current basis orbitals.
///
///          1st argument (vector): indices of orbitals to replace.
///
///          2nd argument (string): name of the data object holding the orbital rotations.
///
///          3rd argument (vector): indices of which rotations to use for the replacement.
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_RotateBasisOrbsSelect_command_function(formic::Archive & arch, formic::AllData & data) {

  if (formic::mpi::rank() == 0 ) {
    formic::of << "   =====================================================================" << std::endl;
    formic::of << "   ======     Replacing Select Orbitals with Chosen Rotations     ======" << std::endl;
    formic::of << "   =====================================================================" << std::endl;
    formic::of << std::endl;
  }

  const int norb = data.userinp.get<int>("nao");

  // read in indices of orbitals that will be replaced
  std::vector<int> repl_vec;
  formic::mpi::read_and_bcast(arch, repl_vec, "failed to read indices of orbitals that will be replaced in formic_RotateBasisOrbsSelect_command_function");

  // read in name of data object storing the matrix of orbital rotation coefficients
  std::string rot_name;
  formic::mpi::read_and_bcast(arch, rot_name, "failed to read name of data object storing the matrix of orbital rotation coefficients in formic_RotateBasisOrbsSelect_command_function");

  // read in indices of the rotations that will be used
  std::vector<int> rot_vec;
  formic::mpi::read_and_bcast(arch, rot_vec, "failed to read indices of the rotations that will be used in formic_RotateBasisOrbsSelect_command_function");

  // call the function that will execute the basis set transformation
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hd, data.wd, repl_vec, rot_vec, rot_name, true);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hc, data.wc, repl_vec, rot_vec, rot_name, true);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_RotateBasisOrbsSelect_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "RotateBasisOrbsAll" command, which rotates the orbital basis.
///
///          1st argument (string): name of the data object holding the orbital rotations.
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_RotateBasisOrbsAll_command_function(formic::Archive & arch, formic::AllData & data) {

  if (formic::mpi::rank() == 0 ) {
    formic::of << "   =================================================" << std::endl;
    formic::of << "   ======     Rotating All Basis Orbitals     ======" << std::endl;
    formic::of << "   =================================================" << std::endl;
    formic::of << std::endl;
  }

  const int norb = data.userinp.get<int>("nao");

  // read in name of data object storing the matrix of orbital rotation coefficients
  std::string rot_name;
  formic::mpi::read_and_bcast(arch, rot_name, "failed to read name of data object storing the matrix of orbital rotation coefficients in formic_RotateBasisOrbsAll_command_function");

  // get indices of orbitals that will be replaced
  std::vector<int> repl_vec;
  for (int i = 0; i < norb; i++)
    repl_vec.at(i) = i;

  // get indices of the rotations that will be used
  std::vector<int> rot_vec;
  for (int i = 0; i < norb; i++)
    rot_vec.at(i) = i;

  // call the function that will execute the basis set transformation
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hd, data.wd, repl_vec, rot_vec, rot_name, true);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hc, data.wc, repl_vec, rot_vec, rot_name, true);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_RotateBasisOrbsAll_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "SetBasisOrbsSelect" command, which sets the coefficients for
///          a selected set of the basis orbitals.
///
///          1st argument (vector): indices of orbitals to replace.
///
///          2nd argument (string): name of the data object holding the coefficient matrix
///                                 for the new orbitals.
///
///          3rd argument (vector): indices of which new orbitals to use
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_SetBasisOrbsSelect_command_function(formic::Archive & arch, formic::AllData & data) {

  if (formic::mpi::rank() == 0 ) {
    formic::of << "   =========================================================" << std::endl;
    formic::of << "   ======     Setting Select Orbital Coefficients     ======" << std::endl;
    formic::of << "   =========================================================" << std::endl;
    formic::of << std::endl;
  }

  const int norb = data.userinp.get<int>("nao");

  // read in indices of orbitals that will be replaced
  std::vector<int> repl_vec;
  formic::mpi::read_and_bcast(arch, repl_vec, "failed to read indices of orbitals that will be replaced in formic_SetBasisOrbsSelect_command_function");

  // read in name of data object storing the matrix of new orbital coefficients
  std::string oc_name;
  formic::mpi::read_and_bcast(arch, oc_name, "failed to read name of data object storing the matrix of new orbital coefficients in formic_SetBasisOrbsSelect_command_function");

  // read in indices of the new orbitals that will be used
  std::vector<int> new_vec;
  formic::mpi::read_and_bcast(arch, new_vec, "failed to read indices of the new orbitals that will be used in formic_SetBasisOrbsSelect_command_function");

  // call the function that will execute the basis set transformation
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hd, data.wd, repl_vec, new_vec, oc_name, false);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hc, data.wc, repl_vec, new_vec, oc_name, false);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_SetBasisOrbsSelect_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "SetBasisOrbsAll" command, which sets the orbital coefficients
///          for the one electron orbital basis.
///
///          1st argument (string): name of the data object holding the coefficient matrix
///                                 for the new orbitals.
///
/// \param[in]     arch         archive to read command arguments from
/// \param[in]     data         data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_SetBasisOrbsAll_command_function(formic::Archive & arch, formic::AllData & data) {

  if (formic::mpi::rank() == 0 ) {
    formic::of << "   ======================================================" << std::endl;
    formic::of << "   ======     Setting All Orbital Coefficients     ======" << std::endl;
    formic::of << "   ======================================================" << std::endl;
    formic::of << std::endl;
  }

  const int norb = data.userinp.get<int>("nao");

  // read in name of data object storing the matrix of new orbital coefficients
  std::string oc_name;
  formic::mpi::read_and_bcast(arch, oc_name, "failed to read name of data object storing the matrix of new orbital coefficients in formic_SetBasisOrbsAll_command_function");

  // get indices of orbitals that will be replaced
  std::vector<int> repl_vec(norb, 0);
  for (int i = 0; i < norb; i++)
    repl_vec.at(i) = i;

  // get indices of the new orbitals that will be used
  std::vector<int> new_vec(norb, 0);
  for (int i = 0; i < norb; i++)
    new_vec.at(i) = i;

  // call the function that will execute the basis set transformation
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hd, data.wd, repl_vec, new_vec, oc_name, false);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_set_or_rotate_basis_helper_function(data.userinp, data.hc, data.wc, repl_vec, new_vec, oc_name, false);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" in formic_SetBasisOrbsAll_command_function.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Add Hamiltonian commands to the formic command system
///
///////////////////////////////////////////////////////////////////////////////////////////////////

void formic::init_hamiltonian_commands() {

  formic::add_command("Psi3Ints", &formic_Psi3Ints_command_function);
  formic::add_command("TextInts", &formic_TextInts_command_function);
  formic::add_command("ReorderOrbitals", &formic_ReorderOrbitals_command_function);
  formic::add_command("WriteTextInts", &formic_WriteTextInts_command_function);
  formic::add_command("ReadBinHam", &formic_ReadBinHam_command_function);

  // add commands for setting the current basis orbitals
  formic::add_command("SetBasisOrbsAll",    &formic_SetBasisOrbsAll_command_function);
  formic::add_command("SetBasisOrbsSelect", &formic_SetBasisOrbsSelect_command_function);

  // add commands for rotating the current basis orbitals
  formic::add_command("RotateBasisOrbsAll",    &formic_RotateBasisOrbsAll_command_function);
  formic::add_command("RotateBasisOrbsSelect", &formic_RotateBasisOrbsSelect_command_function);

  // add command to do an Edmiston-Ruedenberg localization on a given set or orbitals
  formic::add_command("ER_Localize", &formic::ER_Localize_command_function);

  // add commands for Lowdin localization
  formic::add_command("Lowdin_Localize_All", &formic::Lowdin_Localize_All_command_function);
  formic::add_command("Lowdin_Localize_Select", &formic::Lowdin_Localize_Select_command_function);

  // add command to create a new set of integrals with a frozen core
  formic::add_command("Make_FC_Text_Ham", &formic::Make_FC_Text_Ham_command_function);

}
